#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define TRUE    1;
#define FALSE   0;

/* Armazena os elementos para os cálculos do método Jacobi-Richardson*/
typedef struct{
    int J_ORDER;
    int J_ROW_TEST;
    float J_ERROR;
    int J_ITE_MAX;
    float **MA;
    float *MB;
}Jacobi;

/* Armazena o vetor com as soluções encontradas e o número de iterações */
typedef struct{
    float *S_SOLUTION;
    float *S_SOLUTION_AUX;
    int S_ITE;
}Solucao;

/* Faz a leitura dos valores necessários para os cálculos */
void inicializa_Jacobi(Jacobi *J, FILE *file){
    int i, j;

    fscanf(file, "%d", &(J->J_ORDER));
    fscanf(file, "%d", &(J->J_ROW_TEST));
    fscanf(file, "%f", &(J->J_ERROR));
    fscanf(file, "%d", &(J->J_ITE_MAX));

    J->MA = (float**) malloc(J->J_ORDER * sizeof(float*));
    for(i=0; i < J->J_ORDER; i++)
        J->MA[i] = (float*) malloc(J->J_ORDER * sizeof(float));

    J->MB = (float*) malloc(J->J_ORDER * sizeof(float));

    for(i=0; i < J->J_ORDER; i++){
        for(j=0; j < J->J_ORDER; j++){
            fscanf(file, "%f", &(J->MA[i][j]));
        }
    }

    for(i=0; i < J->J_ORDER; i++)
        fscanf(file, "%f", &(J->MB[i]));
}

/* Estritamente Diagonalmente Dominante */
/* Verifica se a soma dos elementos de uma linha, retirando-se a diagonal, será menor que a diagonal */
/* Caso a soma seja igual o maior o método de Jacobi-Richardson irá divergir */
int converge(Jacobi *J){
    int i, j;
    int aux;

    for(i=0; i < J->J_ORDER; i++){
        aux = 0;

        for(j=0; j < J->J_ORDER; j++){
            if(i != j)
                aux += J->MA[i][j];
        }

        if(aux >= J->MA[i][i])
            return FALSE;
    }

    return TRUE;
}

/* Método Jacobi-Richardson */
void Jacobi_Richardson(Jacobi J, Solucao *S){
    int i, j;
    float aux1 = 0.0f, aux2 = 0.0f, erro1 = 0.0f, erro2 = 0.0f;

    S->S_ITE = 0;
    S->S_SOLUTION = (float*) malloc(J.J_ORDER * sizeof(float));
    S->S_SOLUTION_AUX = (float*) malloc(J.J_ORDER * sizeof(float));

	// Cálculo dos resultados
    do{
		// Armazena em um vetor auxiliar o atual vetor de soluções
		S->S_SOLUTION_AUX[i] = S->S_SOLUTION[i];
		
        // Iteração de cada linha
		for(i=0; i < J.J_ORDER; i++){
            S->S_SOLUTION[i] = J.MB[i];

            for(j=0; j < J.J_ORDER; j++){
                if(i != j)
                    S->S_SOLUTION[i] -= S->S_SOLUTION_AUX[j] * J.MA[i][j];
            }

            S->S_SOLUTION[i] /= J.MA[i][i];
        }

		// Cálculo dos erros para critério de parada		
        erro1 = 0.0f;
        erro2 = 0.0f;
        for(i=0; i < J.J_ORDER; i++){
            aux1 = S->S_SOLUTION[i] - S->S_SOLUTION_AUX[i];
            aux2 = S->S_SOLUTION[i];

            if(aux1 < 0)
                aux1 *= -1;
            if(aux2 < 0)
                aux2 *= -1;

            if(aux1 > erro1)
                erro1 = aux1;
            if(aux2 > erro2)
                erro2 = aux2;
        }

        S->S_ITE++;
		
	// Critério de parada é baseado em caso o erro calculado torne-se menor que o erro requisitado,
	// e caso o número de iterações alcance o número de iterações máximo
    }while((erro1/erro2 > J.J_ERROR) && (S->S_ITE < J.J_ITE_MAX));
}

/* Imprimi o vetor de soluções segundo o padrão requisitado */
void imprimi_Solucao(Jacobi J, Solucao S){
    int i;
    float result = 0.0f;

    for(i=0; i < J.J_ORDER; i++)
        result += S.S_SOLUTION[i] * J.MA[J.J_ROW_TEST][i];

	printf("-------------------------------------------------------\n");
    printf("Iterations : %d\nRowTest: %d => [%f] ?= %f\n", S.S_ITE, J.J_ROW_TEST, result, J.MB[J.J_ROW_TEST]);
	printf("-------------------------------------------------------\n");
}

int main(void){
    Jacobi J;
    Solucao S;
    FILE *data;

    data = fopen("data-t3/matriz250.txt","r");

    inicializa_Jacobi(&J, data);

    fclose(data);

	/* 
	 * Caso o sistema seja convergente, é aplicado o método Jacobi-Richardson
	 * caso contrário o sistema é divergente.
	*/
    if(converge(&J)){
        Jacobi_Richardson(J, &S);
        imprimi_Solucao(J, S);
    }
    else
        printf("Sistema linear divergente!");

    return 0;
}

